Descargar datos de internet con tidyquant

Data download with tidyquant

Ana Luisa Espinoza López

Utilizaremos la función tq_get()

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(fpp3)

Registered S3 method overwritten by 'tsibble':
  method               from 
  as_tibble.grouped_df dplyr
── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──
✔ tsibble     1.1.6     ✔ feasts      0.4.1
✔ tsibbledata 0.4.1     ✔ fable       0.4.1
── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
✖ lubridate::date()    masks base::date()
✖ dplyr::filter()      masks stats::filter()
✖ tsibble::intersect() masks base::intersect()
✖ tsibble::interval()  masks lubridate::interval()
✖ dplyr::lag()         masks stats::lag()
✖ tsibble::setdiff()   masks base::setdiff()
✖ tsibble::union()     masks base::union()

library(tidyquant)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
── Attaching core tidyquant packages ─────────────────────── tidyquant 1.0.11 ──
✔ PerformanceAnalytics 2.0.8      ✔ TTR                  0.24.4
✔ quantmod             0.4.28     ✔ xts                  0.14.1── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
✖ zoo::as.Date()                 masks base::as.Date()
✖ zoo::as.Date.numeric()         masks base::as.Date.numeric()
✖ dplyr::filter()                masks stats::filter()
✖ xts::first()                   masks dplyr::first()
✖ zoo::index()                   masks tsibble::index()
✖ tsibble::interval()            masks lubridate::interval()
✖ dplyr::lag()                   masks stats::lag()
✖ xts::last()                    masks dplyr::last()
✖ PerformanceAnalytics::legend() masks graphics::legend()
✖ quantmod::summary()            masks base::summary()
✖ tidyquant::VAR()               masks fable::VAR()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching package: 'tidyquant'

The following object is masked from 'package:fable':

    VAR

library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

library(ggplot2)

You can add options to executable code like this

tq_get_options()

 [1] "stock.prices"       "stock.prices.japan" "dividends"         
 [4] "splits"             "economic.data"      "quandl"            
 [7] "quandl.datatable"   "tiingo"             "tiingo.iex"        
[10] "tiingo.crypto"      "alphavantager"      "alphavantage"      
[13] "rblpapi"           

exports <- tq_get(
1  x = "EXPMX",
2  get = "economic.data",
3  from = "1985-01-01",
4  to = "2025-07-01"
)

exports 

1

El símbolo del dato económico que queremos descargar

2

El tipo de dato y fuente de donde queremos descargarlos, en este caso, datos económicos del FRED

3

Fecha de inicio de los datos aaaa-mm-dd

4

Fecha final

# A tibble: 487 × 3
   symbol date       price
   <chr>  <date>     <dbl>
 1 EXPMX  1985-01-01 1135.
 2 EXPMX  1985-02-01 1117.
 3 EXPMX  1985-03-01 1261.
 4 EXPMX  1985-04-01 1237.
 5 EXPMX  1985-05-01  863.
 6 EXPMX  1985-06-01 1377.
 7 EXPMX  1985-07-01  820.
 8 EXPMX  1985-08-01 1406.
 9 EXPMX  1985-09-01 1016.
10 EXPMX  1985-10-01 1171.
# ℹ 477 more rows

La `tibble` resultante tenemos que convertirla a tsibble y asegurarnos que la variable de fecha tiene el formato correcto

exports <- exports |> 
  mutate(date = yearmonth(date)) |> 
  as_tsibble(index = date)

exports

# A tsibble: 487 x 3 [1M]
   symbol     date price
   <chr>     <mth> <dbl>
 1 EXPMX  1985 Jan 1135.
 2 EXPMX  1985 Feb 1117.
 3 EXPMX  1985 Mar 1261.
 4 EXPMX  1985 Apr 1237.
 5 EXPMX  1985 May  863.
 6 EXPMX  1985 Jun 1377.
 7 EXPMX  1985 Jul  820.
 8 EXPMX  1985 Aug 1406.
 9 EXPMX  1985 Sep 1016.
10 EXPMX  1985 Oct 1171.
# ℹ 477 more rows

1. Convertiremos la variable al formato `yearmonth`
2. Convertimos la tibble a tsibble, indicando que la variable de la fecha es date.

EDA

exports |> 
  autoplot(price)

Ejercicio

1. Realizar un pronóstico a dos años para la serie
2. Entregar un informe en Quarto con el análisis y el pronóstico
3. La métrica de error a utilizar para calcular la precisión del pronóstico será el MAE
4. Entrega: 19-sep

Separar Train

exp_train <- exports |> 
  filter_index(. ~ "2021 Jul")

exp_train

# A tsibble: 439 x 3 [1M]
   symbol     date price
   <chr>     <mth> <dbl>
 1 EXPMX  1985 Jan 1135.
 2 EXPMX  1985 Feb 1117.
 3 EXPMX  1985 Mar 1261.
 4 EXPMX  1985 Apr 1237.
 5 EXPMX  1985 May  863.
 6 EXPMX  1985 Jun 1377.
 7 EXPMX  1985 Jul  820.
 8 EXPMX  1985 Aug 1406.
 9 EXPMX  1985 Sep 1016.
10 EXPMX  1985 Oct 1171.
# ℹ 429 more rows

exp_train |> 
  autoplot(price)

exp_train |> 
  model(stl = STL(price, robust = TRUE)) |> 
  components() |> 
  autoplot() 

exp_train |> 
  autoplot(log(price))

Ahora definimos el test:

exp_test <- exports |> 
  filter_index("2021 Aug" ~ "2023 Jul")

Y el validate:

exp_trainfull <- exports |> 
  filter_index(. ~ "2023 Jul")

exp_val <- exports |> 
  filter_index("2023 Aug" ~ "2025 Jul")

Entrenamiento de modelos

Se calculará el lambda óptimo para box cox con el feature de Guerrero para buscar estabilizar la varianza.

exp_lambda <- exp_train |> 
  features(price, features= guerrero) |> 
  pull()

exp_lambda

[1] 0.02049431

exp_lambda

[1] 0.02049431

En la siguiente sección, se definiran los potenciales modelos para la serie. Se considera que al notar diversos outliers y varianza creciente se utilizará el trend y season multiplicativo.

modelos <- exp_train |> 
  model(
    STL = decomposition_model(
      STL(box_cox(price, exp_lambda) ~ season(window = "periodic")), 
      ETS(season_adjust)
    ),
    Drift = RW(box_cox(price, exp_lambda) ~ drift()),
    ETS = ETS(box_cox(price, exp_lambda)),
    Holt = ETS(price ~ trend("M") + season("M")),
    Damped = ETS(price ~ trend("M", phi = 0.85) + season("M")),
    HW = ETS(price ~ error("M") + trend("M") + season("M")),
    exp_dcmp = decomposition_model(
      STL(box_cox(price, exp_lambda) ~ season(window = "periodic"), robust = TRUE),
      RW(season_adjust ~ drift()),
      SNAIVE(season_year)
    ),
    stlf = decomposition_model(
      STL(log(price) ~ season(window = "periodic"), robust=TRUE), 
      RW(season_adjust ~ drift())
    )
  )

Forecast 2 años

# Forecast 24 meses
fc <- modelos |> 
  forecast(h = "2 years")

# Comparar con el test usando MAE
fc |> 
  accuracy(exp_test) |> 
  arrange(MAE)

# A tibble: 8 × 10
  .model   .type      ME  RMSE   MAE    MPE  MAPE  MASE RMSSE  ACF1
  <chr>    <chr>   <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
1 ETS      Test    922.  1991. 1489.  3.06   5.40   NaN   NaN 0.630
2 Holt     Test    540.  1952. 1537.  1.66   5.62   NaN   NaN 0.604
3 HW       Test    500.  1951. 1550.  1.51   5.67   NaN   NaN 0.608
4 Damped   Test    469.  1959. 1557.  1.39   5.70   NaN   NaN 0.613
5 STL      Test    -11.1 1816. 1565. -0.429  5.83   NaN   NaN 0.654
6 exp_dcmp Test  -1300.  2381. 2021. -5.25   7.72   NaN   NaN 0.692
7 stlf     Test  -1507.  2550. 2140. -6.01   8.18   NaN   NaN 0.707
8 Drift    Test  -1182.  2648. 2151. -4.90   8.22   NaN   NaN 0.521

# Graficar forecast vs test
# Gráfico interactivo usando ggplotly
p <- fc |> 
  autoplot(exp_test, level = NULL) +  # level = NULL quita intervalos de confianza
  autolayer(exp_test, price) +
  labs(title = "Forecast vs Real", y = "Price")

plotly::ggplotly(p)

Se tomó como base el modelo ETS, el cual obtuvo el menor MAE de los modelos propuestos y sobre el se busca añadir elementos al modelo que puedan mejorar el performance del forecast.

ETS_mod <- exp_train |> 
  model(
    ets = ETS(box_cox(price, exp_lambda)~ error("A") + trend("A", , phi = 0.75) + season("A"))
  )

fc_ETS <- ETS_mod |> 
  forecast(h = "2 years") 
  

fc_ETS |> 
  autoplot(exp_test)

fc_ETS |> 
  accuracy(exp_test) |> 
  arrange(MAE)

# A tibble: 1 × 10
  .model .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
  <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ets    Test   477. 1831. 1465.  1.43  5.36   NaN   NaN 0.626

ETS_mod |> 
  gg_tsresiduals()

ets_mod <- exp_trainfull |> 
  model(
    ets = ETS(box_cox(price, exp_lambda)~ error("A") + trend("A", , phi = 0.75) + season("A"))
  )

fc_ets <- ets_mod |> 
  forecast(h = "2 years") 
  
combined <- bind_rows(exp_test, exp_val)
fc_ets |> 
  autoplot(exp_test)

save(fc_ets, file = "fc_AnaLuisa.RData")

Análisis y Conclusión

Analizando las métricas de performance para los distintos modelos planteados para el forecast de la serie, el modelo con menor MAE fue el ETS, con un valor de 1488.526. El análisis de los residuos del modelo ETS, aplicado sobre la serie transformada con Box-Cox para estabilizar la varianza, muestra un buen ajuste a la serie histórica. La mayoría de los residuos oscilan cerca de cero, dentro de un rango aproximado de ±0.2, lo que indica ausencia de sesgo sistemático, y la autocorrelación no revela patrones significativos pendientes de capturar, sugiriendo que los errores se comportan como ruido blanco. La distribución de los residuos es simétrica y asemeja una distribución normal, con valores extremos limitados a eventos anómalos, siendo el más notable el desplome de 2020 por COVID.

En conjunto, esto evidencia que el modelo ETS (“AAA”) captura el movimiento de la tendencia y la estacionalidad a largo plazo, manejando de manera ‘confiable’ la dinámica histórica de la serie, aunque no puede anticipar shocks extraordinarios. Los errores pequeños y homogéneos fuera de los eventos extremos refuerzan su robustez y precisión en condiciones normales.